﻿using OpenCvSharp;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

namespace SSIM
{
    public class SSIMCalculator2
    {
        public double CalculateSSIM(string imagePath1,string imagePath2)
        {
            using var image1 = Cv2.ImRead(imagePath1, ImreadModes.Color);
            using var image2 = Cv2.ImRead(imagePath2, ImreadModes.Color);

            if(image1.Empty() || image2.Empty())
            {
                throw new ArgumentException("无法读取图像文件,请检查路径是否正确");
            }

            if(image1.Size() != image2.Size())
            {
                Cv2.Resize(image2, image2, image1.Size());
            }

            const double C1 = 6.5025, C2 = 58.5225;
            Mat I1 = new Mat(), I2 = new Mat();
            image1.ConvertTo(I1, MatType.CV_32F);
            image2.ConvertTo(I2, MatType.CV_32F);

            Mat I1_2 = I1.Mul(I1);
            Mat I2_2 = I2.Mul(I2);
            Mat I1_I2 = I1.Mul(I2);

            // 计算均值
            Mat mu1 = new Mat(), mu2 = new Mat();
            Cv2.GaussianBlur(I1, mu1, new Size(11, 11), 1.5);
            Cv2.GaussianBlur(I2, mu2, new Size(11, 11), 1.5);

            Mat mu1_2 = mu1.Mul(mu1);
            Mat mu2_2 = mu2.Mul(mu2);
            Mat mu1_mu2 = mu1.Mul(mu2);

            // 计算方差和协方差
            Mat sigma1_2 = new Mat(), sigma2_2 = new Mat(), sigma12 = new Mat();
            Cv2.GaussianBlur(I1_2, sigma1_2, new Size(11, 11), 1.5);
            sigma1_2 = sigma1_2 - mu1_2;

            Cv2.GaussianBlur(I2_2, sigma2_2, new Size(11, 11), 1.5);
            sigma2_2 = sigma2_2 - mu2_2;

            Cv2.GaussianBlur(I1_I2, sigma12, new Size(11, 11), 1.5);
            sigma12 = sigma12 - mu1_mu2;

            // SSIM计算
            Mat t1, t2, t3;

            // t1 = 2 * mu1_mu2 + C1
            t1 = new Mat();
            Cv2.Multiply(mu1_mu2, 2.0, t1);
            Cv2.Add(t1, C1, t1);

            // t2 = 2 * sigma12 + C2
            t2 = new Mat();
            Cv2.Multiply(sigma12, 2.0, t2);
            Cv2.Add(t2, C2, t2);

            // t3 = t1 * t2
            t3 = new Mat();
            Cv2.Multiply(t1, t2, t3);

            // t1 = mu1_2 + mu2_2 + C1
            Cv2.Add(mu1_2, mu2_2, t1);
            Cv2.Add(t1, C1, t1);

            // t2 = sigma1_2 + sigma2_2 + C2
            Cv2.Add(sigma1_2, sigma2_2, t2);
            Cv2.Add(t2, C2, t2);

            // t1 = t1 * t2
            Cv2.Multiply(t1, t2, t1);

            // SSIM map
            Mat ssim_map = new Mat();
            Cv2.Divide(t3, t1, ssim_map);

            // 计算均值SSIM
            Scalar mssim = Cv2.Mean(ssim_map);
            return mssim.Val0; // 灰度图只有一个通道

        }
    }
}
